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Abstract. We present a phase-space method for the Bose-Hubbard model based on the Q- 
function representation. In particular, we consider two model Hamiltonians in the mean- 
field approximation; the first is the standard "one site" model where quantum tunneling is 
approximated entirely using mean-field terms; the second "two site" model explicitly includes 
tunneling between two adjacent sites while treating tunneling with other neighbouring sites 
using the mean-field approximation. The ground state is determined by minimizing the 
classical energy functional subject to quantum mechanical constraints, which take the form 
of uncertainty relations. For each model Hamiltonian we compare the ground state results 
from the Q-function method with the exact numerical solution. The results from the Q- 
function method, which are easy to compute, give a good qualitative description of the main 
features of the Bose-Hubbard model including the superfluid to Mott insulator. We find the 
quantum mechanical constraints dominate the problem and show there are some limitations of 
the method particularly in the weak lattice regime. 



1. Introduction 

The prediction |^ that a Bose-Hubbard Hamihonian could be reahzed with cold Bosonic 
atoms trapped in an optical lattice, and the transition between superfluid and Mott-insulator 
phases observed, led to the development of an experiment in which these predictions were 
verified Q, and has lent urgency to the quest for a computationally simple description of the 
phenomena involved. It is desirable to find a description which can cover both of the two 
regimes of qualitatively different behaviour: 

a) Weak Lattice, when the atoms are delocaUzed, and are thus highly mobile; the system is 
weakly correlated. 

b) Strong Lattice, when the atoms are localized; the system is highly correlated. In this case 
there are two subcases 

i) Commensurate, when the number of atoms per site is integral, and the state is a 
product of number states on each site. In this case there can be no mean-field. 

ii) Incommensurate, when the number of atoms per site is non-integral, and in this case 
there can be a mean-field. 

In the case of an integral average site filling, as the strength of the lattice varies, a transition 
from the delocalized to the localized situation takes place, and one passes from state with 
a non-zero mean-field to a state with no mean-field. Conventionally, if there is no mean- 
field, one speaks of a Mott insulator phase, while if there is a non-vanishing mean-field, the 
terminology superfluid phase is applied. 

The dynamics of the Bose-Hubbard model in states for the strong lattice region presents 
formidable technical difficulties. The so-called Gutzwiller approximation Q has been used in 
several treatments, but because it represents the wavefunction as a product of wavefunctions 
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at different sites, the description of intersite correlations is necessarily very approximate. 
Progress has been made (3] |4l in regions near to the weak lattice regime by using a self- 
consistent Bogoliubov method, in which the statistics of the quantum fluctuations are treated 
by what amounts to a Gaussian ansatz. However, this method cannot be used near the strong 
lattice regime because in this case the statistics are far from Gaussian. 

1.1. Overview 

In this paper we want to introduce approximate methods which should be applicable in both 
the superfluid and the Mott insulator regimes, and to show their efficacy in the very simplest 
method used for the Bose-Hubbard model, that of static mean-field theory. We propose 
a phase space method from quantum optics based on the Q-function representation which 
involves an approximate reparamterization of the Bose-Hubbard Hamiltonian in terms of 
Gaussian variables. The advantage of this methodology lies in the simple description of a 
wide range of quantum states for the system. 

We apply this method to two approximations to the Bose-Hubbard Hamiltonian which we 
denote as the one site and two site models. The one site model treats all intersite correlations 
(ie. tunnelling to adjacent sites) by a mean-field term, which essentially decouples the total 
Hamiltonian into a sum of one site Hamiltonians. This has already been introduced with 
the Hamiltonian given by (|2ji. The two site model extends this formulation by explicitly 
including intersite correlations between two adjacent sites, while treating the interactions with 
neighbouring sites using the mean-field approximation. 

For the one site Q-function approach, we find rather good agreement over all regimes, 
with the advantage that our phase space method can be evaluated very easily — in a matter of 
seconds on any reasonable workstation. The accuracy is usually about 5%, with the qualitative 
behavior being accurately given. The essence of the result is that the ground state is essentially 
that wavefunction which minimizes the classical energy functional subject to the inequality 
constraints given by the uncertainty principle in two different forms. 

These results are compared with an exact numerical solution to the one site model using 
arbitrary one site states. In this case, the solution is found by reparameterizing the energy in 
terms of the number variance and finding the maximum mean-field for which the energy is 
minimized. The formulation is equivalent to the application of the standard Gutzwiller ansatz. 

The Q-function formulation is easily extended to the two site model when lattice 
homogeniety is assumed. The correct description of the two site quantum statistics requires 
the inclusion of additional constraints which are used to determine the ground state solution. 
The results here give a good qualitative description of the overall features of the Bose- 
Hubbard model. However, a shortcoming of the parameterization on two sites is highlighted 
by the failure of the method to correctly predict the Mott insulator phase, as determined by a 
vanishing mean-field, when compared to the results from sections |3lSand|6l 

The two site model is also solved by exact numerical minimization with the assumption 
of lattice symmetry giving a reduced Hilbert space. Here, the ground state results show the 
pertinent features of the Bose-Hubbard model, with a vanishing mean-field throughout the 
Mott insulator phase, which is an improvement over the two site Q-function formulation. 
The results are compared with other reports from the literature, including density matrix 
renormalization group and Quantum Monte-Carlo methods, which yield numerically exact 
results for one dimensional finite lattices. 
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(a) One site model 



(b) Two site model 



Figure 1. Treatment of hopping terms in the Bose-Hubbard model for a lattice of dimension 
d = 1. The dashed arrows indicate those hopping terms tijaiO^^ included by a mean- 
field approximation, whereas the solid arrows indicate those terms included explicitly in the 
Hamiltonian. We consider two formulations: (a) The hopping terms between the ith site and 
it's z = 2d nearest neighbours are treated solely by a mean-field approximation; (b) The 
hopping terms between two adjacent sites i and i + 1 are included explicitly, whereas the 
hopping terms between each of these sites and and their (other) nearest z' = 2d— 1 neighbours 
are included using the mean-field approximation. 



2. Bose-Hubbard Hamiltonian 

In order to set our notation and define terminology, we will summarize the basics of the 
Bose-Hubbard model — first introduced by Fisher et al. |j5l to describe the superfluid-insulator 
transition in liquid ''He, but applicable more generally to any system with interacting bosons 
on a lattice. The simplest form of the model is given when the atoms are loaded adiabatically 
into the lattice, so that they remain in the lowest vibrational state. As we will be formulating 
the model using the Q-function representation, and noting that this distribution is adapted 
to antinormal ordering, we write the Bose-Hubbard Hamiltonian in the antinormally ordered 
form 



where a\ and respectively are the creation and annihilation operators for a boson on the 
ith lattice site. The first term of the Hamiltonian describes the quantum tunelling (hopping) 
between neighboring sites with an amplitude which is only non-zero when i and j are 
adjacent. The second term represents the on-site interaction energy which is taken as repulsive 
with u > 0. We note that this Hamiltonian is easily related to its normally ordered counterpart 
using the commutation relation [a^, a}-] = 6ij. 

2.1. The mean-field approximation 

In the mean-field approximation, one assumes that the effect of the the nearest neighbours on 
the site i is given by a c-number mean-field £ = (a, ) which, by a choice of phase, we take as 
real, and which is independent of i for a homogeneous system. Thus we can write 



where Z is proportional to the nearest neighbor value of Uj multiplied by the number of 
nearest neighbours. As noted in Sachdev's book 151, to solve the system one finds the lowest 
eigenvalue of H„^f on a single site, for which the value of £ must match that of (ui). 

Thus, one assumes a value for the mean-field, computes the ground state eigenfunctions 
of the Hamiltonian, whose mean-field should equal that initially assumed. Unlike the more 
conventional method, we do this at fixed mean occupation of the site, not fixed chemical 




(1) 




(2) 



A phase-space method for the Bose-Hubbard model 



4 



potential, and compute the chemical potential from the energy. Since the solutions of this 
process are well known, it will provide an interesting testbed for the phase space method we 
wish to propose. 




Figure 2. Ground state phase diagram for the Bose-Hubbard model using a self-consistent 
mean-field approach; the mean-field and chemical potential are both shown as a function of 
the relative interaction strength Z/u and mean occupation number n (normal ordered). 



2.2. Self-consistent mean-field results 

The results obtained by this procedure are shown in Fig|2] The usual features of the superfluid 
to Mott insulator transition are present. In particular, in the Mott insulator phase, the mean- 
field vanishes for commensurate occupations when the relative interaction strength {Z/u) is 
below a critical value. There is a corresponding energy gap evident in the chemical potential 
which represents the energy required to add or remove a particle to the system; hence the Mott 
phase is incompressible. When the density is incommensurate or the tunnelling dominates, 
there is a transition to the superfluid phase where the mean-field is non-zero. 

PART I: ONE SITE FORMULATION 



3. Q-functions for one site states 

We first introduce the Q-function parameterization and its application to the one site 
Hamiltonian (|2ji for the Bose-Hubbard model. 



3.1. Phase space methods 

The field of quantum optics leads to various quasi-classical distributions |7l which can be 
used to treat quantum processes using a c-number formalism in terms of coherent states \a). 
For a one-mode system with density operator p, the most widely used of these are given by 

• The Q-function: 

Q{a,a*) = {a\p\a) /t:. (3) 
The Q-function is a quasiprobablity for antinormally ordered operator averages, that is 

(a"(a^)") = y"d2aa"(a*)™g(a,a*) (4) 

It always exists and is positive 
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• The P-function: P(q;, a*) with p = J cPaP{a, a*)\a) {a\, associated with normally 
ordered moments. However, it does not always exist as a positive well behaved function. 

• The Wigner function: W{a,a*), is associated with symmetric operator ordering, and 
always exists, but is not always positive. 

• The positive P-function: P{a, f3), a P-function defined in a doubled phase-space, always 
exists and is positive, but does present some technical difficulties |7|. 

We develop a treatment in terms of the Q-function, since this is always positive and well- 
defined, and is thus a genuine probability density to which we can apply probabilistic 
approximations. It can be used to describe a wide range of states for the Bose-Hubbard model 
which interpolate between the extremes: 

i) Weak lattice: in the case that the hopping dominates (tij very large) the ground state is 
the product of coherent states at each site. 

ii) Strong lattice: in the other extreme, when the hopping is negligible, the ground state is 
the product of eigenstates at each site. If the mean occupation per site is n and [n] denotes 
the integer part of n, then the state at each site is a superposition of number states of the 
form \/A| [n]) + Vl — M + 1)' where A is chosen to give the correct mean occupation. 
These states are not Gaussian. 

The Gaussian or non-Gaussian nature of the statistics in the Bose-Hubbard model is very 
important, and one of the virtues of the Q-function is that it is Gaussian when the quantum 
statistics is also Gaussian Q. 

The Q-functions for the three principal kinds of states are: 

i) Number state: The Q-function for a number state is 



Examples of this distribution appear in Fig|3l (a) and (e). The distribution has a 
characteristic ring shape, which is clearly non-Gaussian. 

ii) Superposition of number states The Q-function for the superposition of number states 
with n — 1 and n atoms in proportions 1 — A : A is 



This is the kind of state expected with an incommensurate filling in the limit of no 
hopping. The Q-function then varies between a ring shaped distribution for A = 0, 1 
to a rather distorted Gaussian distribution at A = 1/2; see Fig|3l(a)-(e). 

iii) Coherent state In the superfluid case with mean filling per site n, where the hopping is 
dominant, we expect an approximately coherent state with parameter /3 — cxp{i6)y/ri, 
for some real 6, and then the Q-function has the Gaussian form 




(6) 




(7) 



This is plotted in Fig|31(f). 
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(a) n = 7, A = 



(b)n = 7, A = 0.15 



(c)n = 7, A = 0.5 




(d) n = 7, A = 0.85 



(e) n = 7, A = 1 



(f) coherent state, /3 = \/6. A : 




Figure 3. Q-function for the superposition of n = 6 and n = 7 states in proportion 1:0, 
0.85:0.15, 0.5:0.5, 0.15:0.85, 0:1; and for a coherent state with (3 = 



3.2. Parametrization of the Q-function 

The probability distributions of Fig|3] can be approximately parametrized by writing the 
random variable a ^ a with 

a={v + (5)e'^ (8) 

where 

a) w is a nonrandom positive real quantity; 

b) (5 is a Gaussian random variable with zero mean and variance a; 

c) is a Gaussian random variable with a mean which is in principle nonzero, but which by 
choice of phase definition can be chosen to be zero, and with variance A. In this case, 
the averages of powers of exp(i6') are given by 

(eP*^) = exp ( VA/2) = c^'. (9) 

Using properties of Gaussian variables this approximation leads to the values of the 
antinormally ordered moments: 

(10) 
(11) 



i^a) = vc 



(aa) = (v^ + a)c^ 



(12) 
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(aaa^) = {v^ + 3(Jv)c 
{aaa^a^) = (u^ + 6v^a + 3(7^) 



(13) 
(14) 



Note that the first of these moments ilOi represents the mean-field. A straightforward 
extension of this this parameterization to L sites can be given by the substitution 
{vk + 5k)e^^'' for k = 1,2, ...,L; we make use of this for the two site problem presented 
in section|5| 

Validity of parameterization: values of the parameters in particular cases To fit the kinds 
of distribution in Fig|3lwe determine the parameters v, a and c by fitting the moments (a), 
(aa^) and {aaa'^a''). These are given in detail in [Appendix A| where we show that for a 
wide variety of states we get very tolerable approximations. Thus we can expect a good 
quaUtative description of the system for an arbitrary lattice strength. However, we note that 
the parameterization is least accurate in the case of an equal superposition of number states 
which we will find reflected in the ground state phase diagram in the strong lattice regime. 

3.3. Determination of the ground state 

Using the Hamiltonian (|2j and the moments J10> - J14> . the average energy in the Q-function 
representation is 



To find the ground state solution at a given interaction strength Z/u and mean occupation n, 
this quantity is minimized with respect to the free parameters v and c, with the bounds w > 
and < c < 1 as permitted by the Q-function parameterization. 

3.3.]. Constraints The functional il5\ is a classical distribution that has no local minima so 
that the minimum must occur on the boundary. However the system should be governed by 
the quantum mechanical nature of the problem, which we reintroduce in the form of two 
constraints on the minimization procedure, one derived from a restriction on the number 
variance, the other from the uncertainty relation for the conjugate variables of phase and 
number 

The most significant of the constraints enforced by the quantum mechanical nature of the 
problem are as follows. 

Fractionality constraint If the mean occupation per site is non-integral, the variance of the 
occupation must be nonzero, and this has a major bearing on the problem, which we shall 
formulate precisely. We will be considering cases where the mean occupation per site is 
fixed, so that 



E{n, c, v) 



(15) 




n > 1 



(16) 



so that 



n = + (T > 1. 



(17) 



Setting 



(18) 



and using ( I10> - J14> the variance of the site occupations is 



= 2n^ -n- 2v^. 



(19) 
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In n is non-integral, with fractional part 5n, then the minimum variance for a given n occurs 
when only the two occupation numbers which bracket the value n are represented, and then 
this gives a variance of 5n{\ — Sn), leading to the constraint 

var[A^] > Sn{l - Sn). (20) 

and using ( I19> we can write this as 



< u"^ < 



— ( n + 5n{l — Sn) 



(21) 



Phase-number uncertainty relations The uncertainty principle enters when we consider 
that phase and number are conjugate variables. We can write a rigorous uncertainty 
relationship using 

X^^; Y=^ (22) 

so that we have the commutation relations 

[N, X]^ - iY; [N, Y] = iX (23) 
from which follow an uncertainty relation 



{SY^)var[N] > ^{X)^ 



(24) 



In our formulation we have chosen a zero mean phase so that (Y) = (see ( IA.13» and thus 
{SY'^) = {y^). Using ( IA.IBMaTTtI . and antinormally ordering all the products involved, we 
find 

^ (25) 



(<5r2) = |(i-c^)-^ 



and thus 



(1 



(2n^~n~2v^) > -wV. 



(26) 



3.3.2. Results Minimizing (I15> with the constraints i2l\ and i26\ leads to the ground state 
phase diagram shown in Fig|5] These results compare well with the exact solution shown in 
FigEl particularly with confirmation of the Mott insulator phase at commensurate occupations 
for a strong lattice (small Z/u), and of the superfluid phase for a weak lattice (large Z/u). 

The results are dominated by the interplay between the constraints, which is made 
evident by considering Fig@] In the strong lattice regime, the ground state solution is given 
at the intersection of the fractionaUty constraint (I21> which determines the upper bound for v, 
and the phase-number uncertainty constraint (I26> . shown by the solid curve. Specifically for 
Z/u ~ 0.5, and for an integer n ~ 1, this occurs at c = and the resulting vanishing mean- 
field is indicative of the Mott insulator phase; for a non-integer ri = 1.5 the solution occurs 
for c 7^ so that there is a non-zero mean-field corresponding to a superfluid component. 
Conversely for a weak lattice, with n = 1.5 and Z/u = 10, the solution is determined solely 
by(|26j. 

The fact that the quantum mechanical constraints dominate the problem is not unexpected 
as the phase transition between the Mott insulator and superfluid is driven by quantum 
fluctuations. 

We note that for the Q-function results at half-integer values of n, the mean-field remains 
constant over some range of Z/u in the strong lattice regime. This contrasts with the 
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(a) n = 1, Zju = 0.5 



(b) n = 1.5, Zju = 0.5 



(c) n = 1.5, Z/ii = 10 



Figure 4. Level curves of energy 1151 in the Q-function representation for three different 
cases are indicated by the dashed curves, with the energy decreasing for darker regions; the 
permissable parameter space for v and c is given by the bounds -u > and < |c| < 1 (for 
purely illustrative purposes we allow c to be negative here because the problem is symmetrical 
about c = 0) where v is further bounded by 1211 . Since 1151 has no local minima, the minimum 
energy occurs within the shown bounded region and subject to the constraint 1261 derived from 
the phase-number uncertainty relation which is shown by the solid curve. The solution in each 
case is indicated by a circle. In the first case (n = 1, Zju = 0.5) this occurs for c = so that 
the mean-field is zero; in the second case (n = 1.5, Zju = 0.5) this occurs for a non-zero 
c so that there is a mean-field. In both these cases, the solution saturates the upper bound for 
V. Conversely, in the third case (n = 1.5, Z/u = 10) where the lattice is weak, the solution 
occurs within the bounded region. 



numerically exact results obtained by a self-consistent mean-field approach as shown in Fig|2 
where the mean-field is a monotonically increasing function of Z/u. This difference arises 
from the weaker approximation of the Q-function parameterization in the case of an equal 
superposition of number states as shown in Fig|3](c), and quantified in [Appendix A| 




Figure 5. Ground state results for Bose-Hubbard model with Q-function representation; the 
mean-field and chemical potential are shown as a function of relative interaction strength Z/u 
and (anti-normal ordered) occupation n. 
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4. Exact numerical minimization of one site problem 

The inequalities \26\ and \2\\ restrict the mean-field to a certain maximum value, but this 
maximum may not be optimal. Let us then pose the question: What is the largest value of the 
mean-field for a given mean number and variance. This value determines the minimum energy 
and permits the following basic strategy for finding the ground state phase diagram. Firstly, 
reparameterize the mean-field in terms of the mean number occupation and number variance; 
and secondly, minimize the energy with respect to the number variance at a fixed relative 
interaction strength zt/u and mean occupation. The problem is then essentially reduced to a 
minimization of energy with respect to a single variable - the number variance. 



4.1. Formulation using one site states 

Here we give an exact numerical procedure for determining the ground state of the Bose- 
Hubbard model in the one site mean-field approximation. A one site quantum state can be 
written as a superposition of number states 

|s)=^c„|n) (27) 

71 

where the usual normalisation condition (1 = {s\s)) applies. In this basis the mean-field is 

{X) = {s\'^^\s) ^^[£VTicl_,Cr] (28) 

The Hamiltonian (|3 leads to the expression for the one site energy in the mean-field 
approximation 

{E) ^ -Z {Xf + u{2 + var[iV] + + 3?t) (29) 

where n = {a' a) is the (normal ordered) mean number and var[Af] = — is the number 
variance of atoms on each site (note the normal and anti-normal ordered means are related by 
(a^a) = (aa^) — 1). Since the chemical potential term has been omitted in ([0, we treat the 
system using the canonical ensemble by taking a fixed number of atoms in the lattice, that is 
by explicitly using the constraint of fixed mean occupation. 

Clearly, the energy functional i29\ is then minimised when the mean-field {X) is a 
maximum and subject to the constraints of normalisation and fixed number mean and variance. 
In terms of the one site state i27\ these constraints are given by 

— ^ |c„p — 1 = (normalization) (30) 

n 

(j)2 = |c„p n — n = (fixed mean) (31) 

n 



El 



' n - = (fixed variance) (32) 



Moreover, since the constraints are independent of phase, the mean-field ( I28> is maximized 
by choosing the set of Cn to be real and positive at fixed |c„|. To deal with the problem 
of constrained optimisation we apply the method of Lagrange multipliers; the mean-field is 
maximized when 

O^^-iA^ + ^^ + ii.^ (33) 
dcn 2 dcn 2 dcn 2 dCn 
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Note the scalar factors that appear before the Lagrange multipliers (A, /i and v) are included 
for convenience. Equation ( I33> along with equations ( I28> and (|30}-(|32} then lead to 

V«c„_i + Vn+l c„+i + (/in + vr?)cn = Ac„ (34) 

This is a Hermitian eigenvalue problem which has only one eigenvector solution with the 
set c„ all positive, since the solutions form an orthonormal basis (solutions with some c„ 
vanishing may cause ambiguities in principle however). 




Figure 6. Mean number occupation as a function of at = —500 for tlie one site model. 
The Lagrange multiplier ^ plays the role of (but is not equal to) the chemical potential; the 
horizontal regions are indicative of the system near the incompressible Mott insulator phase 
where dn/d^ = 0. 



4.2. Numerical Solutions 

Numerical solutions of ( I34> determine the maximum mean-field in terms of var[A^] and n. 
By choosing a fixed Z/u and n the calculation of the ground state energy then amounts to a 
direct minimization of \29\ with respect to the variance var[Af]. 

We outline the procedure as follows. First, for an appropriate range of the Lagrange 
multipliers (/i, v) the eigenvector solutions of ( I34t are used to calculate corresponding values 
of n, var[Af] and the maximum mean-field {X). A fixed value of v determines a n{^) 
curve (see for example Fig|6j which is inverted using linear interpolation - with an additional 
optimization step - to give /i on a uniform set of n points. That is, each value of n determines 
a curve in ji and v from which the dependence of {X) on var[iV] is determined from solutions 
to ( I34t . Since this relationship is determined on a finite set of points, it is necessary in practice 
to use cubic interpolation to determine {X) for arbitrary var[iV]. 

Calculations of {X) as a function of var[Af] are shown in FigEl for commensurate 
and incommensurate mean site occupations; the corresponding curves which arise from the 
uncertainty relation M6\ are shown to be very similar. Note that for commensurate mean 
occupations the mean-field and variance simultaneously approach zero, corresponding to the 
Mott insulator phase (a pure number state). In contrast, in the incommensurate case neither 
the mean-field nor variance approach zero; even in the strong lattice regime there is a non- 
zero superfluid component, corresponding to a superposition of two number states. Using 
these curves, the energy \29\ can be minimized with respect to var[iV]. This leads to ground 
state results that are the same numerically as those of the self-consistent mean-field approach 
shown in Fig|2] 
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Although the difference in the two sets of curves shown in Fig appears slight, it is 
entirely responsible for the notable difference between the exact one site calculation (FigO 
and the approximate phase-space calculation (Fig|5jl- In the incommensurate case, the fact 
that the exact curves drop rapidly as var[Af] 1/4 is responsible for the sloping behaviour 
of the mean-field arches shown in Fig|2]as opposed to the level behaviour shown in Fig|5] 



Exact 

Commutation relation 

3.5 1 . . . . . 1 3.5 




var[AG var[Ar| 



Figure 7. Mean-field versus variance for commensurate (a) and incommensurate (b) mean site 
occupations for the one site case. 

Numerical calculations were performed using the MATLAB software package. The state 
space was truncated with n < 40 as this provided a good trade-off between computational 
efficiency and accuracy (ie. this truncation introduced negligible error to solutions in the 
region of interest n < 10). Note, to cover a sufficient n and var[iV] range, {n, ly) were chosen 
from the triangular region with —800 < < —0.01 and < /i < —40 v. The minimization 
procedure used the MAIL AB fininbnd (bounded minimization) function, with a lower bound 
for the number variance determined by the fractionality constraint i20\ . 

4.3. Calculating the phase transition boundary 

To check the validity of the one site formulation near the transition point, it is useful to 
calculate the the position of the phase boundary between the superfluid and Mott insulator 
phases. This can be done approximately by a perturbation method where the ground state 
energy is determined by states near the transition point. We consider a normalised state of 
commensurate occupation n near the transition point 

\s) = VT^ln) +y/X{e''^^n-l) +€'"'111 + 1)) (35) 

where A is a real and small; A = corresponds to the Mott insulating phase with a 
commensurate occupation n. States with A > characterize the superfluid phase where 
the mean-field is non-zero. The corresponding one site energy j29> is given by 

E= -A(l-2A)(Vncos6'i + VWTTcos6i2)^ 

+ u{2 + 2A + + 3n) (36) 

Note that in scaling by the factor 1/zt this form of the energy is dimensionless and we have 
defined the ratio u = u/zt for the relative interaction strength. The ground state is determined 
by the minimum energy with respect to variations in the free parameters: A, 6i and 02. Clearly 
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this occurs for 9i ^ 62 = (the result 61 ~ 62 ~ gives the same ground state energy and 
is equivalent to a change in the sign of the mean-field) and when 

— = = -(1 - 4A)(\/Wcos6'i + V^lTTcosez)^ + 2u (37) 

The transition to a Mott insulator phase occurs in the limit A ^ when the mean-field goes 
to zero. Applying this condition, the critial point for the transition is then given by 

uc = 5(V^ + VWTT)2 (38) 

For values of the relative interaction strength above the critical value u > Uc the ground state 
solution saturates the bound at A = and the system remains in the Mott insulator phase. 

Comparison to other work In order to compare our results with those reported elsewhere we 
define the parameter 

Uc = 2{u/zt)c (39) 

to describe the relative interaction strength at the critical point. The factor of 2 is included 
because the on-site interaction term u in the Hamiltonian Q is equivalent to the term U/2 
which has been used elsewhere (see |4l and |8l for example). 

Equation ( I38t then becomes Uc = (v^+ + 1)^, an expression also found elsewhere 
Ell^l; this yields the following values for the transition: Uc — 5.83 for n ~ 1, Uc = 9.90 for 
n = 2,Uc = 13.93 for ?T = 3. 

PART II: TWO SITE FORMULATION 

The one site formulation necessarily neglects correlations between neighbouring sites 
in the mean-field approximation. We can improve on this by using a Hamiltonian that 
explicity includes hopping between two adjacent sites while still treating interactions with 
other neighbouring sites with the mean-field approximation. In particular consider two 
adjacent sites labelled by 1 and 2 corresponding to sites i and i + 1 respectively as shown 
in Fig. n^b). Following the one site case Q we can write the Bose-Hubbard Hamiltonian in a 
two site mean-field approximation as 

-fftwosite = -t aial + a2a\ + {2d-l) (^^£{ai + a\) + ^£{a2 + al)^ 

+ u {aiaia\a\ + a2a2a\a\) (40) 

for a homogenous lattice of dimension d where £ = {ak) (for k 7^ 1, 2) is the mean-field 
representing the interaction with nearest neighbours, which as with the one site model, by a 
choice of phase is taken as real. 

5. Q-function parameterization on two sites 

In this case the Q-function parameterization becomes 

ai ^ ai = {vi + Si)e'^' (41) 

a2^a2 = {v2 + hV^^ (42) 
We define the intersite correlation operators 

6 = aia|; &^ = a2a| (43) 

The corresponding moments are 

(6) = (fot) = (i;it.2+u.)5 (44) 
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where we have set 

w = {6162) (45) 

(46) 

and assumed that the phase and number correlations are independent; although this 
assumption is possibly not strictly valid, without it the formalism becomes significantly more 
complicated. 

The two site Hamiltonian (I40> leads to the expression for the average energy 

E = (-ff two site) 

= - t [2{viV2 + w)g + {2d - l){vlcl + wfc^)] 

+ -inf - 2vf + in\ - 2v\ (47) 

which has been normalised by dividing through by u and defining the relative interaction 
strength t = t/u. 

5.1. Constraints 

Following the one site case, to find the ground state solution, the energy ( I47> is minimized 
subject to constraints arising from any applicable uncertainty relations. In particular, in 
addition to the constraints (I21> and \2b\ that have already been introduced for the one 
site statistics, we can derive further constraints that account for the two site statistics by 
considering commutation relations between bilinear operators. 

Variances and covariance In deriving the necessary constraints for the two site problem, we 
first consider the following operators 

A^i + N2 

N = (48) 

Ni - No 

M = — (49) 

2 

The corresponding variances for each of these operators are given by 

var[7V] = l/4(var[A^i] + var[iV2] + 2cov(iVi, iVs)) (50) 
var[M] = l/4(var[iVi] + var[iV2] - 2 cov(iVi, TVa)) (51) 
where the number covariance function is given by 

cov(7Vi,iV2) = {N1N2) - {Ni){N2} 

^ 2{2viV2W + w^) (52) 
with the one site variances, var[Afi] and var[Af2], given analogously to jl9t . 

Commutation relations We define the following operators 

b + b^ 

(53) 

W = — — (54) 
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{V) is the intersite correlation and gives a measure of particle exchange, or tunnelling, 
between adjacent sites. Eqs. ( I53l l and i54l have the following commutation relations 

[V, M] = iW (55) 

[M, W] = iV (56) 

[W, V] = iM (57) 

Quantum mechanical constraints Using ( I55> -( I57> and noting that (W) — 0, we can write 
the corresponding uncertainty relations 

((^2) - {Vf) var[Af] > \{[SV, SM]+f (58) 
var[M]{W^) > ^(vf + \{[5M,5W]+f (59) 
_ ^vf) (W) > \{Mf + \{[5V,5W] + f (60) 



with 



= ^^Kn2+cov(7Vi,7V2))-^^ (61) 
W = (wi 1-2+^)5 (62) 
(Vt/2) = + cov(7Vi, iV2)) - ^^^^ (63) 



Applying homogeneity These constraints further simplify when homogeneity is assumed. In 
particular, if we assume the ground state is symmetric to the 1^2 interchange, we can set 

vi ~ V2 ~ V (64) 

ci = C2 = c (65) 

ni = 722 ~ n (66) 

Clearly then 

(A/) = (67) 
Moreover, under the above interchange the following moments are equal 

(6tiV2) = {hNi) (68) 

{hN2) = (fe^iVi) (69) 
It follows that the anti-commutators in equations (I58>-(I60> are given by 

[5V,5M]+ = Q (70) 

[5M,5W]+ = Q (71) 

[5V,5W]+^Q (72) 

Two site phase-number uncertainty relation Similarly, the commutation relation (i = 1, 2) 

[N,Y,]^\iX, (73) 
leads to the constraint 

var[7V](F,2) > (74) 
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zVu ° ^ n zt/u -\ 



(a) Mean-field, d=l (b) Intersite correlation, d=l 





(c) Mean-field, d=2 



(d) Intersite Correlation, d=2 




(e) Mean-field, d=3 (f) Intersite Correlation, d=3 

Figure 8. Ground state phase diagram for two site Bose-Hubbard model using the Q-function 
representation in the mean-field approximation. The mean-field and intersite correlations are 
shown as a function of the mean (anti-normal ordered) occupation n and the relative interaction 
strength zt/u. 
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Two site fractionality constraint A lower limit on the two site variance var[iV] can also be 
given by considering two site states with minimum variance. Such a state has the form 

\s) = ^/Z^|[n], [n]) + ^^(1 [n], [n] + 1) + |[n] + 1, [n])) 

+ V7lW + l,W + l) (75) 

where [n] = n — Sn is the integer part of n, 6n being the fractional part. We then have 

var[7V,;] = (76) 



2 V 2^ 

cov(7Vi,iV2) = -/3V4 (77) 
so that we can write the minimum variance 

var[7V]„„„ = - /3) (78) 

Using ( 15 Ot and applying homogeneity this leads to the fractionality constraint for the two site 
variance 

- in - + 2v^w + > \(3{1 - (3) (79) 
We note the following two cases 

(i) < i5n < 0.5 : In this case 7 = and [3 = 2Sn. 

(ii) 0.5 < 6n < 1 : In this case a = and /3 = 2(1 - 6n). 

5.2. The ground state solution 

To find the ground state solution we assume homogeneity and minimize the two site energy 
J47> at a given occupation n and relative interaction strength t with respect to the free 
parameters v, c, w and g. This is done subject to the constraints (12 It . \26\ . ( l58> - (l60t . MA\ 
and ( I79> . The parametrization requires that 

w > (80) 

Additionally, by noting the variables c and g both take the form (e"'^), the following bounds 
apply 

0<c<l (81) 
< g < 1 (82) 

5.3. Results 

The results of this calculation are shown in FiglS] The solutions reproduce the general features 
of the Bose-Hubbard model in the superfluid limit and for incommensurate occupations in 
the limit of zero tunneling. However, in the case of a commensurate mean occupation, the 
Q-function parametrization on two sites does not correctly demonstrate the existence of the 
Mott insulator phase. Specifically we find that for a three dimensional lattice, the mean-field 
only vanishes in the limit of zt/u ^ {]. Moreover, in the one and two dimensional cases, the 
predicted Mott insulator phase only occurs for very small values of zt/u. See sections|3]|4]and 
|6lfor comparison; in particular, the results from section|6lshow a well defined Mott phase with 
a vanishing mean-field. Improved agreement can be seen between the intersite correlations 
calculated by the Q-function representation and the exact numerical results of section|6l In 
both cases the intersite correlations vanish in the limit of zero tunneling for commensurate 
occupations. 
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6. Exact numerical miniinization of two site problem 

Following the one site formulation discussed in section lTTI we apply a two site generalization 
of the exact numerical calculation using arbitrary two site states. The formulation uses the 
assumption that the lattice is translationally invariant so that we can apply symmetry for the 
two sites in the ground state. 

The Hamiltonian j4m leads to an expression for the two site energy 

(£^)twosite = - 2t [{V) + {2d-l){Xf] 

+ 2u (2 + var[7V] + + 3 n) (83) 
where V is the intersite correlation operator given by ( I53> and 

l = i(a, + a|) (84) 

is the mean-field operator, which by homogeneity, does not depend on the site index. 

Note that in the limit that {V) — ^ {X)"^ the two site energy tends to twice the one site 
energy. This corresponds to the case where the two site state is factorizable IV'12) iV'i) 1^2) 
and the one site behaviour is returned. 

In general, for n and m atoms on sites 1 and 2 respectively, the two site quantum state 
can be written as 



= E CN.,p{\\{N+p),\{N-p)) + \\{N~p),\{N+p))) 

+ Cn,o\IN,^N) (85) 

N,p=0 

where we have introduced the quantum numbers N = n + m and p = \n — m\ < N to 
account for the lattice symmetry c„ m = Cm.n (we should note this symmetry argument can 
only be appUed when the ground state exhibits homogeneity). The correlation function is 

{V) = (a.a]) 



Cn, mCn+1, m-l\/m{n + 1) 

n. m 

J2 Cn.,pCn.p+2^{N~p){N+p+2) 
\5p^iY.^Ct,,if{N + l) 



N. p>0 



(86) 



N 

and the mean-field is 



^ C'w.P [Cn+Up+i\/N+p+2 + Cn+i.p-iVN-P+2 



V2 



N,p>0 



^-5^,0 5Z^^iv,oCAf+i,iViV+2 (87) 
V 2 „ 
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From ( 1^ it is evident that, for a given n, d and t/u, the ground state solution can be found 
by finding the maximum kinetic (hopping) energy 

F = {V) + {2d-l){X,f (88) 

subject to the constraints 

(i) normalisation, 1 = (s|s) 

(ii) fixed mean occupation (per site), n ~ ^{a\ai + a\a2) 



(iii) fixed variance (per site), v? = ^{a[aia{ai + a2a20'2'^2) 



'i\aia\ai + a\a2a\t 
In terms of the symmetric state ( I85> these constraints are respectively given by 

2 E '^w,p + E *^w,o -1 = (normalisation) (89) 

N,p>0 



N 



^N,pN +\Y. '^n.qN \-n = (fixed mean) 

N,p>0 J 

I E ci,iN-'+p') + lY.cUN']-^^o 

N,p>0 N J 



(90) 



(fixed variance) 

Using the method of Lagrange multipliers, F is then maximized when 

1 . ^(|)^ d4>2 d(j)3 







dF 



-X- 



OCn.p 2 OCn.p dCN.p 



dCM: 



(91) 
(92) 



6.0.7. Vectorizing the problem We can write this set of simultaneous differential equations 
in matrix form when the set of coefficients Cat, p is expressed as an ordered vector To truncate 
the problem space in terms of a maximum two site occupation A^max (with a corresponding 
Pmax), it is convenient to use an ordering for the quantum numbers p and N with increasing p 
(for each N value) within increasing N. That is, 

C = [C'o,0: Cl,li ^2,0, C2^2, C34, 6*3^3, Cifi, Ci^2, C4_4 . . . C N,-n,^^,p„^.^^] (93) 

The correlation and mean-field can then be written respectively in terms of the vector 
quadratic forms 

(V^) = C^AvC (94) 

{X) = C^AxC (95) 

where the construction of matrices Ay and Ax follows directly from equations ( 15 3> and ( I84> . 
In particular, 

\ 

































1 
































2V2 
























































2 


2^3 



























































2v/6 































4 






















































(96) 



A phase-space method for the Bose-Hubbard model 



20 




Figure 9. A contour plot of the mean occupation n with respect to the Lagrange multiphers 
H and u reveals the Mott lobes for commensurate occupations, shown for n = 1. 2, 3. where 
the solutions to 11041 are degenerate with respect to ^, indicating the incompressibility of the 
phase. The Mott insulator phase does not exist for incommensurate occupations as shown by 
the remaining curves at n = 0.5, 1.5, 2.5 and 3.5 from left to right in the figure. The dashed 
lines represent possible (non-degenerate) trajectories for the Mott insulator phase as discussed 
in section l6. 1.11 
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(97) 



01 = C^5C -1 = 

02 = SDC-n = Q 
(j)3 = lC'^ SD'C-^=0 

where S, D and D' are diagonal matrices with elements given by 

Snn'pp' ~ Snn' Spp' (2 — Spo) 
Dnn'pp' = Snn' Spp' N 



D' 



NN'pp' 



(98) 
(99) 
(100) 

(101) 
(102) 
(103) 
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Figure 10. Mean-tield as a function of number variance for commensurate (a) and 
incommensurate (b) occupations. In contrast to the one site approximation, in the 
commensurate case the Mott insulator phase with vanishing mean-field can occur for a non- 
zero number vaiiance. 



with the usual Kronecker delta function. Using these definitions, equation (|92} can then be 
written in matrix form as 



-1/2 .0-1/2 



where 



A- 
B 



'AS- 

-■ Av H 



Au 



{2d-l){X){A 



X 



A^ 



X) 



(104) 

(105) 
(106) 
(107) 



This can be made linear in u (and therefore C) when the mean-field term {X) in ( I106> is 
treated as a free parameter x; the resulting eigenvalue equation is solved iteratively in x so that 
it is self-consistent with the calculated mean-field {X) = C^AxC. The eigenvector solution 
u which maximizes F at each iteration is nonnegative; in fact the nonnegative eigenvector, 
the maximizing solution, belongs to the maximum eigenvalue as we now show. 



6.0.2. Determining the maximising solution To see why this is the case we need to invoke 
the Perron-Frobenius theorem for nonnegative matrices 1101 . We first note that the r x r matrix 
G ~ S'~^/^AS'^^/^ + S'^/^i?S'~^/^ in equation ( I104> is nonnegative every where except along 
the leading diagonal where negative values are possible when either ji or: v (or both) are 
negative. It is then always possible to construct a matrix G' = G + jlr that is nonnegative 
everywhere for a sufficiently large 7, being the r x r identity matrix. 

The eigenvectors of G" and G are clearly the same with the eigenvalues of G" given 
by A'; = Ai + 7 where A; is the ith eigenvalue of G corresponding to the zth eigenvector. 
The Perron-Frobenius theorem states that the spectral radius p of a nonnegative matrix is an 
eigenvalue corresponding to a nonnegative eigenvector^. That is, p{G') is the eigenvalue with 

f The spectral radius for an n X n matrix A with eigenvalues A; (1 > i > n) is defined as p{A) = max| Ai|. 
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(a) Mean-field, d=l (b) Intersite correlation, d=l 




(c) Mean-field, d=2 (d) Intersite correlation, d=2 




(e) Mean-field, d=3 (f) Intersite congelation, d=3 

Figure 11. Ground state results for the Bose-Hubbard model in the two site mean-field 
approximation. The mean-field and intersite coiTelations are shown as a function of the mean 
site occupation n and the relative interaction strength zt/u. A non-zero mean-field represents 
the superfluid regime. The mean-field is zero below a critical value (zt/u)c for commensurate 
occupations, indicating the onset of the Mott insulator phase. The intersite correlations are 
non-zero everywhere except for commensurate occupations in the limit of no tunnelling where 
zt/u = 0. The results are qualitatively similar for 1, 2, and 3 dimensional lattices. However, 
a feature of the d = 1 phase diagram is the weak suppression of the mean-field at half-integer 
occupations, which is seen more seen more clearly in Fig |13l This may be an artefact of the 
two site approximation and requires further study. 
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(a) Mean-field, d = I 



(b) Mean-field, d = 2 



(c) Mean-field, d = 3 





(d) Correlation, d =1 



(e) Correlation, d = 2 



(f) Correlation, d = 3 
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n=2 
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zt/u 




(g) Variance, d = 1 



(h) Variance, d = 2 



(i) Variance, d = 3 



Figure 12. d = 1, 2, 3 results: mean-field, intersite correlation and variance for commensurate 
occupations of n = 1, 2, 3. The dashed vertical lines indicate the corresponding positions 
(from left to light) of the phase boundary between the superfluid and Mott insulator phases as 
calculated using the pertubation theory outhned in section lS!3l 



the desired nonnegative solution. We can relate this to the maximum eigenvalue of G when 7 
simultaneously satisfies the following two conditions; 

(i) 7 > |min(G)| so that G" is nonnegative 

(ii) 7 > |min(Ai) I so that all X[ > and max(Ai) + 7 is necessarily the spectral radius of G' 

For any bounded matrix G it is always possible to select such a value of 7; we can therefore 
associate the correct solution for Gn.p with the maximum eigenvalue of ( I104> which returns 
a self-consistent mean-field. 
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6.1. The ground state solution 

The general procedure for finding the ground state solution follows the method used in the one 
site formulation (see section|3. In particular, at a given mean occupation n, ( I104l i determines 
a locus of solutions which can be expressed as a function of number variance. Therefore, at 
a given t/u and n, the two site energy can be reparameterized in terms of number variance; 
finding the ground state phase diagram then only requires the minimization of the energy 
with respect to variance. In practice to perform this procedure efficiently, it is useful to 
first calculate a set of solution states covering a suitable variance range, which act as initial 
conditions in the final minimisation step. 

In the results presented here, this was achieved by fixing n and initially choosing a large 
value for the variance corresponding to the superfluid regime (specifically var[iV] = n for 
n < 1 and var[Af] =n/2 for n > 1). The corresponding state was found by solving equation 
(fT04l using the initial values xq ~ 5, /^to = 2 and i^q = —1, which was found to be sufficient 
for the range of parameters considered here. A lower variance was then selected and its 
corresponding state was calculated, by using the previous solution to determine the initial 
conditions. By iteration it was then possible to find solutions efficiently over a broad variance 
range. 

However, although this general procedure works well for most of the parameter space, 
it fails to converge (to a sufficient accuracy) in two limiting cases: for Mott-insulator states 
where the mean-field is zero and there is degeneracy; and in the limit of the lowest possible 
variance, var[iV] = 5ri{l — Sn), corresponding to the localized region where the Lagrange 
multipliers diverge at the solution. It is therefore necessary to consider these cases separately. 

6.1.1. Calculating Mott states For commensurate occupations and below a critical value for 
the variance, the system is in the Mott insulator state corresponding to a vanishing mean- 
field. In this case, the above approach is numerically unstable for two reasons. Firstly, the 
algorithm fails to correctly converge to a zero mean-field since the search direction cannot be 
determined. Secondly, the system states become highly degenerate in ji in the Mott insulator 
state (corresponding to a incompressible phase) which further prevents convergence (with 
respect to fj). 

These problems are circumvented respectively by the following adjustments to the 
general procedure. Firstly, explicitly setting the free-parameter x to zero in the equation 
| |104> . ensures that any solution is self-consistent in the mean-field; this can be seen to be 
case by considering how the mean-field term appears in the eigenvalue equation (see ( I94> and 
J106> in particular). Secondly, setting ji = —r v for a suitably chosen quantity r, restricts the 
parameter space to the degenerate region, so that n remains fixed while var[7V] ^ in the 
limit V ^ oo. These trajectories are illustrated by the dashed curves in Fig|9] By enforcing 
this relationship between ^ and v, we can therefore probe the Mott insulator states accurately. 

6.1.2. Solutions in the limit of zero tunneling Considering ( I83> in the limit of zero tunneling 
where t/u 0, the ground state is clearly determined by states of the lowest possible 
variance. In this case, we found the method of Lagrange multipliers leads to poor convergence 
as the Lagrange multipliers diverge at the solution. We therefore consider a more direct 
method whereby we minimize the two site energy in terms of a reduced parameter space, 
which correspond to the lowest allowed variance (var[7V] = 5n{l — 5n)) for a given mean 
occupation n. 

Recall the lowest variance state is given by equation M5V . in this basis, the hopping 
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(a) (b) 



Figure 13. The mean-field (a) and intersite correlation (b) using a reduced basis of three states 
in the limit of zero tunelling (zt/u = 0) for cZ = 1, 3. The methodology is outlined in section 
16.1.21 A feature of the d = 1 case is the weak supression of the mean-field around half-integer 
occupations. The d = 2 results, which been omitted for clarity, are very close to the d = 3 
case. 

energy ( I88> is given by 

F = P/2{[n] + 1)(1 + {2d - 1)(V^ + V7)') (108) 

The constraints (normalization, fixed n and fixed variance) can be given as a set of linear 
equations 

a + /3 + 7 =1 
+ ([n] + i)/3+ ([n] + l)7 =n 
[rT]^a+ ([n]V[n] + i)/3+ ([n]^ + 2[n] + l)7 = v&i-[N]+n^ (109) 

which when written as Ax = b, admits solutions of the form x = xq + k nullspace[A] where 
we have taken 





a 




[n\ + 1 — n 




Xo = 


P 







(110) 




7 




n — [n] 





as a solution of the homogenous system Ax ~ 0. There is an additional requirement on k 
that the coefficients are bounded with 0<a<l,0</3<2 and < 7 < 1. For a given 
occupation n and dimension d, we can then maximise the hopping energy ( I108t with respect 
to the remaining free parameter k. The ground state phase diagram is then calculated easily 
for the case where zt/u = 0; the resulting mean-field and intersite correlation are shown in 
figures ll^T a) and lOf b) respectively. 

6.1.3. Numerical Solutions Numerical calculations were performed using a state space 
truncated with N < 20, corresponding to 121 coefficients Cm.p- The MATLAB function 
fsolve was used to find the self-consistent solution, for a given n and var[Af], by supplying the 
function with suitable initial values xq, fiQ and i/q for the mean-field and Lagrange multipliers. 
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6.2. Results 

Using the calculated states with maximum tunneling energy, the two site energy was then 
minimized with respect to var[A^] for a given n and Z/u. The results of this calculation are 
shown in Figs. II If al-Cf) for arbitrary mean occupations. To investigate the superfluid to Mott- 
insulator transition more clearly, these results have also been reproduced for commensurate 
occupations in Figures [T^ a)-(i). 

This ground state phase diagram is qualitatively similar to the results for the one site 
model (Fig|2}; the most important feature of both models is the vanishing of the mean- 
field for commensurate occupations above the critical transition point Uc, corresponding to 
the Mott insulator phase. Note that the Mott insulator to superfluid transition is clearly 
reproduced for commensurate mean occupations, an improvement over the corresponding 
two site Q-function approximation (section |5} which can be attributed to the exact treatment 
of the intersite correlations in the present treatment. We now consider a simple perturbative 
expansion from which the position of the phase transition boundary can be easily calculated 
for the two site formulation. 



6.3. Calculating the phase transition boundary 



Following the one-site case, we consider the following states near the transition with 
commensurate occupation n: 

lot , 



\s) = yJl-a-2P |n,n) + 



(|n- 1) + |n + l,n - 1)) 



2 I I"'" 



1) + |n- l,n) + 1) + \n+l,n) 



(111) 



This state is a superposition of the minimum set of number states that give a non-zero mean- 
field and intersite correlation. The coefficients have been chosen to satisfy normalization 
and a fixed commensurate occupation n. The corresponding variance on one site is given by 
var[iV] = a + /3. Note that \s) represents a perturbation in a and j3 from the Mott insulator 
phase. With this state the intersite correlation ( I86> is given by 



and the mean-field (l84l is 



/3, 



2/3)a ^/n{n + l) + |(2n + 1) 



(112) 



(113) 



In the ground state, the phase transition occurs for the minimum two site energy ( I83t : clearly 
this occurs with respect to the state ( II 1 1> when dExijda = and dEi2/dj3 = are both 
satisfied. Moreover the phase transition occurs for /3 = where the mean-field is zero. That 
is, 



= 



= 



dE 






da 


13=0 


(1) 


dE 






9/3 


13=0 


(7) 



v/n(l + n)(l-2a) 
v/2a(l - a) 

- (fj-i(l + 2n) 

i(2d - + VWTT)2(\/i~l^ 

y^2an{n + 1) 



(114) 



a 



(115) 
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d 


n = l 


n = 2 


n = 3 


1 


4.75 


8.02 


11.26 


2 


5.56 


9.43 


13.27 


3 


5.71 


9.69 


13.63 


10 


5.82 


9.88 


13.90 


10^ 


5.83 


9.90 


13.93 



Table 1. Critical value U c calculated for two site mean-field approximation using a 
perturbation expansion around the Mott Insulator state 



Numerically solving these equations leads to the results shown in tabled Note that in the 
non-physical limit of d — > oo the variance at the transtion is var[Af] = a = 0. Considering 
equation M \5\ the transition then occurs at Uc = 2{u/2dt)c = (V^ + + 1)^, the same 
value as in the one-site model. This is not unexpected as the intersite correlation {V) becomes 
negligible compared with the mean-field contribution in the Hamiltonian (I40> . 

There are some notable differences between the one and two site mean-field 
approximations, which we discuss here. In particular, the inclusion of intersite correlations 
in the two site Hamiltonian, results in a non-zero number variance even in the Mott insulator 
phase where the mean-field vanishes; this can be seen in Figs. [TOlandE] This contrasts with 
the one site results (see Fig. where the mean-field and number variance are identically 
zero. Moreover, for the two site model the calculated transition points Uc, as shown in table 
[2 occur at lower values than for the one site model. In particular, for n = 1 and d = 1, the 
two site model yields Uc ~ 4.75, whereas the one site model gives Uc = 5.83. 

These differences are most pronounced for the c? = 1 case where the mean-field 
approximation is no longer valid, as has also been noted elsewhere iSj . At higher dimensions 
(d ~ 2, 3), where the mean-field approximation is more accurate, the number variance is close 
to zero in the Mott phase and the calculated transition points for the one and two site models 
are close. It is expected that in the (non-physical) limit of infinite lattice dimensionality, the 
one and two site results will converge. 

6.4. Comparison to other work 

The results for the two site Hamiltonian given in the previous section clearly demonstrate a 
partial shift to lower values U c (ie- a weaker lattice) when compared to the one site results. 
This shift can be attributed to the inclusion of intersite correlations in the model. To compare 
with our results with those reported elsewhere it is useful to refer to table|2] 

For the d ~ 1 case, numerically exact schemes, with Density Matrix RenormaUzation 
Group (DMRG) and Quantum Monte-Carlo (QMC) simulations, predict a much lower value 
for the phase transition point Uc in the Bose-Hubbard model. Using the DMRG technique 
with the infinite-system algorithm, the transition point for n = 1 and d = 1 has been reported 
as Uc = 1-68 1 1 I ! and Uc = 1-81 1121 . Using the more accurate finite-system algorithm, 
the same transition has been reported as U c = 1-92 |TJj and U c = 1-68 1141 . Similarly, for 
QMC simulations without lattice disorder, the transition has been reported as Uc = 2.33 1151 . 
Also, pertubative expansions using defect states have yielded similar results to the that of the 
QMC technique in one dimension. In particular, the above transition has been reported as 
Uc = 2.33 kl61l andUc = 2.04 El- 
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Method 


Position of phase transition bound- 
ary for d = \, n = 1 (unless other- 
wise stated) 


Authors 


Reference 


Mean-field approximation 
using second order pertur- 
bation expansion 


C7c = 2n + 1 + ^{2n + 1)2 - 1, 
n^l:Uc = 5.83 


Oosten et al. 




Mean-field theory using 
variational approach 


Uc = 2.0 


Amico and Penna 


118) 


Mean-field theory using 
variational approach 


n = 1: [/c = 5.83 


Sheshadri et al. 


El 


Perturbative expansion us- 
ing defect (particle/hole) 
states 


d = l,n = I: Uc = 2.33; d = 2, 
n = 1: C/c = 3.68 


Freericks and Monien 


m 


Perturbative expansion us- 
ing defect (particle/hole) 
states 


T T C\ r\ A / ill* 

U c ~ 2.04 (see table m reference 
for further calculated values at dif- 
terent a and n) 


Freericks and Monien 


(ni 


Quantum Monte-Carlo 


Uc = 2.33 (without disorder) 


Batrouni and Scalettar 




DMRG with infinite- 
system algorithm 


Uc = 1.68 


Pai et al. 


lu 


DMRG with infinite- 
system algorithm 


Uc = 1.81 (without nearest- 
neighbour interactions) 


Kuhner and Monien 


m 


DMRG with finite-system 
algorithm 


Uc = 1.68 (without nearest- 
neighbour interactions) 


Kuhner et al. 


m 


DMRG with finite-system 
algorithm 


Uc ~ 1.92 (for zero disorder) 


Rapsch et al. 


|T3l 



Table 2. The calculated boundary of the Superfluid to Mott-insulator phase boundary as 
determined from the Bose-Hubbard model as solved using various methods in the literature. 



7. Conclusions 

We have presented two model Hamiltonians for the Bose-Hubbard model using two different 
mean-field approximations. The simpler one site model treats all intersite correlations using a 
mean-field approximation that decouples the problem into the sum of one site Hamiltonians. 
The second model is a two site extension where intersite correlations between two adjacent 
sites are explicitly included while treating interactions with neighbouring sites using the 
mean-field approximation. Each model has been tackled using two methodologies: a 
treatment in terms of a Q-function representation; and a numerically exact method using either 
the one or two site states. 

In the case of the one site Hamiltonian, we find the Q-function representation agrees well 
with the numerically exact treatment, but can be solved at a fraction of the computational cost. 
For the two site Hamiltonian, the Q-function gives a good qualitative description but does not 
give a clear Mott insulator to superfluid transition due to limitations of the parameterization 
on two sites. It is encouraging, however, that the Q-function approach gives the intersite 
correlations accurately when compared with the two site exact solution. In contrast to the 
Q-function approach, the two site exact solution yields a well-defined phase transition for 
commensurate mean occupations. In this case, the critical relative interaction strength U c is 
smaller (corresponding to a weaker lattice strength) for all commensurate values of n than for 
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the one site model; this shift is most pronounced for a one dimensional lattice. 

What emerges from the Q-function approach is that the quantum mechanical constraints 
play a primary role in determining the ground state results for the Bose-Hubbard model. This 
is most clearly illustrated by the one site formulation where the results are determined by the 
uncertainty relation in two forms: a restriction on the minimum variance permitted by the 
quantum state; and a relation between the number and phase fluctuations. In particular, the 
onset of the Mott insulator phase is characterized by supression of number fluctuations at the 
lower bound of the variance constraint. 

Exact numerical results with QMC or DMRG simulations on finite one dimensional 
lattices by other authors indicate the actual position of the superfluid to Mott insulator 
transition occurs for an even lower value of U c than is calculated by our two site model. 
It is expected that if we were to extend our formalism to explicitly include the interactions 
between three or more sites, the predicted value Uc would shift to lower values in line with 
these other treatments. However, this remains a computationally difficult problem due to the 
dramatic increase in the size of the Hilbert space of the problem when more sites are included. 
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Appendix A. Values of parameters in particlar cases. 

To fit the kinds of distribution in Fig|3lwe determine the parameters v, a and c by fitting the 
moments (a), (aa^) and {aaa^a^). 

Parameters for a number state 

The exact values for the quantities ( I10H14> in the case of a number state |m) are 








m + 1 



(A.l) 
(A.2) 
(A.3) 
(A.4) 
(A.5) 




leading to the choice 



V 



a 




(A.6) 



(A.7) 
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~ 1/4 as m — > oo 
c =0. 

The two quantities ( IA.3IA.4t are given exactly by this choice. 

Parameters for a superposition of two number states 
For an equal superposition of the kind Q that 



(A.8) 
(A.9) 



n + 



The exact computation of the quantities ( I10H14> yields in this case 



[aa) 
(aa) 



= 



i(?i + l)V^ 



(aaa) a)) 



{n + lf 

The first, second and last of these are fitted exactly, leading to the approximate values 

(aa) w + a)c^ 
1 



(u^ + 3fTw)c 



16(n2 + n/2- 1/8 

- 2v/n2 + n/2- 1/8 



(A. 10) 
(A.ll) 
(A. 12) 

(A. 13) 

(A. 14) 
(A.15) 
(A. 16) 
(A. 17) 

(A.18) 



(A. 19) 



Even for rt = 1 these are very tolerable approximations; ( IA.18> is diminished by comparison 
with (aa^) by a factor of « 0.045, compared to the exact value of 0. For large n \A.\9\ 
becomes equal to the true value JA.16> . and even for n = 1 the true and the model results 
differ by less than 8%. 

Parameters for a coherent state 

For a coherent state \ j3), where /3 = m is taken to be real, the quantities take the form 

(a) = 13 (A.20) 

1 (A.21) 

(A.22) 



aa' 
aa) 



aaatfl^) = /J'' + 4/3^ + 2 



(A.23) 
(A.24) 
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from which we can compute 



„4 = ^ /32 + (A.25) 



c =-~l-7^ as/3^cx3 (A.27) 

It is easy to check that for large [3 the quantities ( IA.22IA.23l are correctly given, and they are 
of course correct for (3 — 0. For /3 = 1, we find that (aa) = 0.8 instead of the exact value of 
1, and (aaa^) = 2.84, instead of the exact value of 3. 



Appendix B. A rigorous lower bound on the energy 

As a spinoff from our methodology of section|3l we find that we can also develop a rigorous 
lower bound for the energy as a function of the mean-field, and we can evaluate the predictions 
given by assuming the lower bound is equal to the ground-state energy. 

To see this we note the mean energy, which is given by (I29> in the mean-field 
approximation, is parametrized by number variance var[Af] and the mean-field {X) . Then we 
know that var[A^] satisfies the two inequalities ( I20> and (I24> . which we can combine together 
as 

var[iV] > max - ^n)^ ■ (B-l) 

We can develop a lower bound on (F^) by noting that 

(r2)=„_i_(X2)<r^-i-(X)^ (B.2) 
so that we can deduce from iB.H that 

var[iV] > max (- ^l!—— , - 6n)) , (B.3) 

and that the mean energy i29i satisfies the lower bound 

{E)> - Z{Xf + u(2 + n^ + 3n + max( i^^L——,Sn{l-Sn)]\. (B.4) 

I J) 

This bound is parametrized by the mean-field {X) only, so the rigorous lower bound is given 
by the value, X,„i,i(u, n), of {X) which gives the minimum value, £'i„in(w, n), the right hand 
side of iBA\ . and subject to the restriction that {X) < ^/n — 1 (the equality holding for a 
coherent state). These are easy to evaluate, and resulting chemical potential and mean-field 
are plotted in Fig lBll The results are less accurate than the phase space method (see Fig|5j, 
but are still surprisingly good. 
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Figure Bl. Ground state phase diagram calculated using rigorous lower bound on energy 
in Q-function representation 



